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We report the critical point for site percolation for the "explosive" type for 2D square lattices 
using Monte Carlo simulations and compare it to the classical well known percolation. We use 
similar algorithms as have been recently reported for bond percolation and networks. We calculate 
the "explosive" site percolation threshold as pc = 0.695 and we find evidence that "explosive" site 
percolation surprisingly may belong to a different universality class than bond percolation on lattices, 
providing that the transitions (a) are continuous and (b) obey the conventional finite size scaling 
forms. We do not attempt to determine the order of the explosive transition. Finally, we study and 
compare the direct and reverse processes, showing that while the reverse process is different from 
the direct process for finite size systems, the two cases become equivalent in the thermodynamic 
limit of large L. 



I.INTRODUCTION 



Random percolation (RP) [iMZl is a well studied model 
for phase transitions in statistical physics and condensed 
matter physics. It consists in randomly occupying sites 
(site percolation) or bonds (bond percolation) on a given 
lattice with probability p. Neighboring occupied sites 
merge to form clusters. As p increases, new clusters are 
formed or coalesce, until at a critical value, Pc^ a giant 
component emerges, transversing the system, the well- 
known infinite percolating cluster. 

Classical percolation and all of its variants (continuum 
percolation, site-bond percolation, bootstrap percola- 
tion, invasion percolation, etc [l|), are known to be con- 
tinuous transitions, exhibiting a divergent correlation 
length ^ at Pc- Thus, there was a stupendous surprise 
that recently, Achlioptas et al [sl , proposed a new method 
for the occupation of sites which produces "explosive" 
transitions: when filling sequentially an empty lattice 
with occupied sites, instead of randomly occupying a site 
or bond, we choose two candidates and investigate which 
one of them leads to the smaller clustering. The one that 
does this is kept as a new occupied site on the lattice 
while the second one is discarded. This procedure con- 
siderably slows down the emergence of the giant compo- 
nent, which is now formed abruptly, thus the term "ex- 
plosive" . 

Although it is not clear if "explosive" percolation is a 
first-order transition |9 Mll| o r continuous transition with 
a very small (3 exponent [12| , it is certainly a very sharp 
transition with a lot of interesting and unusual proper- 
ties and has stimulated considerable interest. Follow- 



ing these first announcements new work emerged |l2Lll3| 
that doubted the result of the Achlioptas problem being 
a first-order phase transition, thus making this problem a 
hot contentious issue. Recently, applications of the "ex- 
plosive" percolation transition for "real" systems have 
been proposed [13. [l5| . Up to now the systems investi- 
gated include networks and lattice bond percolation. 
In this paper we extend these investigations to include 



now site explosive percolation and its finite-size hystere- 
sis properties. We find that site explosive percolation 
indeed exhibits a very sharp transition with /3 exponent 
"effectively" zero {(3/iy ^ 0.001) for the case of sum rule 
and c:^ 0.04 for the case of product rule. We investigate 
two variants of a reverse "Achlioptas process (API and 
AP2)" and find that the direct and the reverse processes 
exhibit a hysteresis loop, which surprisingly vanishes in 
the themodynamic limit. 

The rest of the paper is organized as follows. In Sec. II 
we describe the site percolation algorithm for the direct 
and the reverse processes. Moreover, we describe the 
tools we used to examine the transition and the behavior 
of the system. In Sec. Ill we present our results and 
discuss their significance. Finally, in Sec. IV we sum up 
our main conclusions. 



II.METHODS 

In the present paper we use Monte Carlo simulations 
for the site percolation problem on L x L square lattice 
with periodic boundary conditions. The algorithm for 
the case of sum rule (APSR) proceeds as follows: 

(1) Initially, we start from an empty lattice. We ran- 
domly occupy one single site. 

(2) Next, we randomly select a trial unoccupied site, say 
A. 

(3) We calculate the size sa of the resulting cluster in 
which A belongs to. 

(4) We remove the trial unoccupied site A and randomly 
select a trial unoccupied site B, different from A. 

(5) We calculate the size 5b of the resulting cluster in 
which B belongs to. 

(6) In case sa < sb-, site A is permanently occupied and 
site B is discarded. In case sb < sa^ site B is perma- 
nently occupied and site A is discarded. In case sa = sb^ 
we randomly select and permanently occupy either A or 
B discarding the other. Each time, the number of occu- 
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pied sites t is incremented by one. 

(7) We repeat steps (2)- (6) until the entire lattice is cov- 
ered. For each "time step" we monitor the size of the 
largest cluster Smax- 

In Fig [T] we present an illustrative example of the algo- 
rithmic process, as has been previously studied. Next, 
we implemented the product rule (APPR). In this case, 
we calculate the product of the sizes of the different clus- 
ters which will merge after the placement of the new site, 
A or B. The only changes in the above algorithm are at 
steps (3) and (5). In order to compare the "explosive" 
percolation transition with the classical site percolation 
transition, we have simulated random site percolation, 
using the efficient Newman-Ziff algorithm ||^]. 
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FIG. 1: Achlioptas Process according to the sum rule(APSR) 
for site percolation . White cells correspond to unoccupied 
sites while colored cells correspond to occupied sites. Dif- 
ferent colors (red, green, gray,blue) indicate different clusters, 
(a) We randomly select two trial unoccupied sites (yellow), 
noted by A and B, one at a time. We evaluate the size of the 
clusters that are formed and contain sites A and B, sa and 
SB respectively. In this example sa — 10 and sb — 14. (b) 
According to AP, we keep site A which leads to the smaller 
cluster and discard site B. 

Subsequently, we simulated two variants of reverse "ex- 
plosive" percolation processes, in order to study possible 
hysteresis phenomena associated with "explosive perco- 
lation". We have implemented the following algorithm 
to simulate the reverse process API: 

(1) Initially, we start from a fully occupied lattice. 

(2) Next, we randomly select two trial occupied sites, say 
A and B. We remove them temporarily from the lattice. 

(3) We place again one of the previous sites, say A and 
calculate the size sa of the resulting cluster in which A 
belongs to. 

(4) We go back to step (2) and repeat (3) for the case of 
site B, calculating the size of the resulting cluster in 
which B belongs to. 

(5) In case sa > sb^ site A is reoccupied and site B 
is permanently discarded. In case sb > sa^ site B is 
reoccupied and site A is permanently discarded. In case 
5a = ^5^, we randomly select and reoccupy either A or B, 
permanently discarding the other. Each time, the num- 
ber of occupied sites t is reduced by one. 

(6) We repeat steps (2)- (5) until the entire lattice is 
empty. For each "time step" t, we monitor the size of 



the largest cluster Smax- 

In the above algorithm we did not consider periodic 
boundary conditions. In Fig [21 we present an illustrative 
example of the algorithm used. The details are given in 
the figure caption. We start with a fully occupied lattice. 
The resulting clusters are shown in Fig |2(d)[ 
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FIG. 2: Reverse Achlioptas Process (API) for site percolation 
according to the sum rule. Blue is for the occupied sites while 
white for the unoccupied sites. Initially, the lattice is fully 
occupied. (a) An instance of the process. We randomly choose 
two trial sites (yellow), noted as A and B, and remove them 
from the lattice, (b) The clusters formed after the removal, 
(c) We place site A again in the lattice and calculate the size 
of the cluster in which it belongs, sa = 16. (d) We do the 
same as before for the case of site B and calculate sb = 26. 
We remove site A which leads to the formation of the smaller 
cluster and keep site B. 

The reverse process AP2 is a slight modification of the 
above algorithm for API. At step 5 we remove the site 
which maximizes the resulting cluster. Thus, for AP2 
step 5 should be: (5 -AP2) In case sa < sb^ site A is 
reoccupied and site B is permanently discarded. In case 
Sb < sa^ site B is reoccupied and site A is permanently 
discarded. In case sa = sb, we randomly select and re- 
occupy either A or B, permanently discarding the other. 
Each time, the number of occupied sites t is reduced by 
one. The rest of the algorithmic steps remain identical 
to API described above. 

Note that the main difference in the selection rule ap- 
plied in the direct and in the reverse API algorithm is 
that in the direct case, we add the site which minimizes 
the resulting cluster, while in the reverse API case, we 
remove the site which minimizes the resulting cluster. 
Also, for clarity, we performed the opposite procedure in 
the reverse case (named reverse AP2 in Fig 7(b)), i.e. we 
removed the sites which maximized the resulting cluster. 

To track the behavior of the system, we have stud- 
ied the order parameter Pmax defined as the ratio of the 
sites that belong to the largest cluster Smax to the total 
number of sites in the lattice. For continuous transitions. 
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Pmax follows the Scaling relation [J, ^ : 

Pma. = L-''^''F[ip-p,)L'^''] (1) 

in the vicinity of Pc. For the present analysis, we assume 
that the order parameter of the transition follows the 
conventional well-known finite size scaling. Recently, it 
has been suggested [ill that "explosive" transitions are 
continuous, but with an unconventional finite size scal- 
ing. We do not examine this possibility in the present 
manuscript. Moreover, we monitor the standard devia- 
tion of the size of the largest cluster Smax^ defined as 



\/ < S^^^ > < Smax >^ 



(2) 



11, 16| as a criterion to 



which has recently been proposed 
discern between continuous and "explosive" transitions. 

We also employed the A criterion proposed in [H which 
quantifies the transition width of a system with size N . 
It is defined as A = ti — to, where is the lowest value 
for which Smax > 0.5N and to the lowest value for which 
Smax > a/A^. Achlioptas et al. stated that for continuous 
(i.e. random) percolation, this quantity should scale as 
A/N for AP and A/N const for RP as TV ^ oo. 
However, this simple scaling relation does not seem to 
hold for bond percolation on square lattices 

The data given in the following figures are averages 
over 1000 realizations. 



III.RESULTS AND DISCUSSION 

In this section we present the results obtained 
from the calculations. In Figs 3(a)[ |4(a) and |5(a) 
we plot the order parameter Pmax as a function of 
p, for the case of RP, APSR and APPR models, 
respectively, for five different system sizes, namely 
L = 200{black),400{red),600{green),800{blue) and 
1000{wine). For the case of RP model, we have re- 
covered the correct critical exponents {P/i^ — 0.106 and 
l/u 0.75). A comparison between Figs 3(a)|4(a) and 
5(a) shows that Pmax changes more abruptly in the case 
of AP models, as expected from an "explosive" transi- 
tion. Moreover, in the APSR case, we observe that at 
p :^ 0.695 all lines for the plotted system sizes seem to 



be crossing roughly at the same point. This p value coin- 
cides with our estimations of the transition point Pc for 
the APSR model, using the method of the second largest 
cluster as well as the value p at which the x per lattice site 
(x/A/") reaches its maximum value. Using standard finite 
size scaling techniques, we find that P/u 0.001. More 
precisely, Eq. [Dimplies that L^^'^Pmax = F[{p-pc)L^^''] 
and, thus, curves of L^^^ Pmax versus p for different sys- 
tem sizes should cross at a single point dX p = Pc- We 
have used our Monte Carlo data for Pmax vs p for 5 differ- 
ent sizes L, to create 5 functions gi = L^^^Pmax and used 
minimization software to determine the values j3/v and 



p which minimize the sum of the squares of the pairwise 
differences K = K{p/v,p) = Y.._^. {gi - gjf . 

Ideally, K{l3/v^Pc) should be exactly zero, but since 
numerical errors cannot be avoided, we consider that the 
values oi j3/v and p that minimize A provide an accu- 
rate estimation of the actual values oi j3/v and p^ Such 
a small value for the P/v as the one observed for the 
APSR, is very unusual for a continuous transition. Of 
course, it cannot be determined whether a transition is 
actually discontinuous or continuous with an extremely 
small /3, as already discussed in ref [lH. It is, however, 
very interesting that our calculated value fS/v is different 
from the values calculated by Ziff [l6| and Radicchi et al. 
[loj for bond percolation. They also could not reach a 
sharp conclusion regarding the nature of the transition. 



In contrast, in [12|, ll3| numerical evidence is provided 



for the continuity of "explosive" percolation transition. 
Our findings indicate that in case the "explosive" tran- 
sition is continuous, the site "explosive" percolation and 
bond "explosive" percolation belong to different univer- 
sality classes. This result is both intriguing and unex- 
pected, since in the classical case, it is well-known that 
the critical exponents do not depend on the structure of 
the lattice (e.g., square or triangular) or on the type of 
percolation (site, bond or even continuum) [2|. We must 
emphasize however, that Grassberger et al. |13| suggest 
that "explosive" percolation follows an unconventional fi- 
nite size scaling. 

In Figs 3(b)||4(b)] and |5(b)[ we plot the normalized stan- 
dard deviation x/N oi the size of the largest cluster Smax^ 
for the same five system sizes L as above, for the RP, 
APSR and APPR case, respectively. In all cases, x/^ 
exhibits a sharp maximum at p^ However, in the RP 
case, we observe that the maximum value is a decreasing 
function of L, while in the cases of APSR and APPR, 
it is nearly constant for large systems (see also Fig |6(b)| 
below for APSR). 

As discussed above. Fig [5] shows our results from the 
product rule simulations (APPR). Using the minimiza- 
tion technique previously described, we calculated the 
value of Pc to be Pc = 0.756 ± 0.006 in accordance to [l^ 
within error bounds. As for j3/v we calculated it to be 
p/u = 0.04 ±0.02. 

In Fig 6(a)[ we plot the quantity A/N as a function 
of the system size L, both for RP (red dots) and APSR 
(black squares) models. We observe a power-law scal- 
ing for both cases, with an exponent ujrp = —0.27 and 
CjOap = —0.39. As expected, and in agreement with pre- 
vious results reported by Achlioptas et al. for the case of 
random networks lH and by Ziff for bond percolation on 
lattices [sl, we find that the exponent for the AP case is 
much smaller than the exponent for the RP case, indicat- 
ing a sharper transition for the former. Again, we observe 
that the exponent uap for the site "explosive" percola- 
tion is different from the exponent for the case of bond 
"explosive" percolation reported by Ziff in jof, providing 
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FIG. 3: Random Percolation (RP).(a) Plot of 
function of p for different linear sizes L. (b) Plot of x/^ as a 
function of p for the same sizes L as in (a). It is evident that 
the peaks of the curves decay as L increases. Colors are: L = 
200 (black), L = 400 (red), L = 600 (green), L = 800 (blue), 
L lOOO(wine) 



an additional indication that site and bond "explosive" 
percolation belong to different universality classes. 

Next, we study how the quantity x/^ scales with the 
system size L for the classical and the "explosive" case. 
For the RP case, we expect, by a simple scaling argu- 
ment, that the quantity x/^ decays as L'^l^^~^. It is 
known that the variance of the order parameter Praax is 
related to the susceptibility exponent 7 and scales at Pc 
as vaT[Pmax) - ^^'''^ Thus, vav^Smax) - The 
standard deviation of the quantity Smax^ X-> is the square 
root of the variance and, thus, x ^ L^l^^^^ . And, finally, 
x/A/" ~ IjII'^^-^ ^ For the RP process, the values of 7 and 
V are known exactly and they are 7 = 43/18 and = 4/3. 
This leads to the conclusion that x/A^ ~ — 5/48 —0.1. 
In Fig |6(b)[ we plot x{Pc)/N as a function of L for RP 
(red dots) and for APSR (black squares). Points are 
Monte Carlo simulation results, while the dotted line is 
the theoretically expected slope for the RP case. We did 
observe that the RP case data show a power-law decay, 
while for the APSR case, we observe an initial increase 
oi x/N followed by a plateau for large system sizes. 
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FIG. 4: Achlioptas Process with sum rule (APSR). (a)Plot 
of Pmax as a function of p for different linear sizes L for site 
percolation. We observe that all plotted lines pass through a 
single point at p 0.695. Mor eover, we see a more abrupt 

(b)Plot of x/Aas; 



jump of Pmax compared to|3(a)l (b)Plot of x/A as a function 
of p for the same sizes L as in (a) . Peaks are almost constant 
Colors are: L = 200(black), L = 400(red), L = 600(green 
L = 800(blue), L = lOOO(wine) 



Finally, in Fig [71 we present the results for the direct 
and the reverse AP processes (API and AP2) described 
in Sec. II. In Fig |7(a)[ we plot 

Pmax ^S a function of 
p for a 700 x 700 lattice, both for the direct and the 
reverse API case. We observe that for a finite system, 
there exists a hysteresis loop. Of course, no such loop 
exists in the classical RP case, where the process is fully 
reversible. We are interested in determining the scaling 
properties of this loop. We should emphasize that in 
contrast to thermal processes, we cannot associate the 
observed hysteresis loops with metastability and first or- 
der transitions. We merely use the "hysteresis" concept 
as a means to investigate and compare different but re- 
lated processes. In Fig |7(b)[ we plot the area Ai^ sur- 
rounded by the curves corresponding to the the direct 
and reverse processes for various sizes L as a function 
of 1/L. Black squares are for the reverse API and red 
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FIG. 5: Achlioptas Process with product rule (APPR). (a) 
Plot of Pmax as a function of p for different linear sizes L for 
site percolation, (b) Plot of x/A^ as a function of p for the same 
sizes L as in (a). Colors are: L = 200(black), L — 400(red), 
L = 600(green), L = 800(blue), L = lOOO(wine) 
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dots are for the reverse AP2, as described in Sec. II. 
The fitting suggests that AH - (l/Lj^-^^s 
mer and AH - (l/L)^-^^^ for the latter. Thus, we find 
that, although in finite systems, the reverse processes are 
different from the direct process, they become equivalent 
in the thermodynamic limit (L oo). 

IV.CONCLUSIONS 



FIG. 6: (a) A /AT as a function of L for RP (red dots) and 
APSR (black squares). Both scale as a power-law, with L~°'^^ 
for RP case and L'^'^^ for APSR case, (b) Plot of x{Pc) /N as 
a function of L for RP (red dots) and APSR (black squares). 
Blue dashed line is for the theoretical prediction of the vari- 
ation of x{Pc)/N for the case of RP. We observe a power-law 
decay for the RP case. In the APSR case, x{Pc)/N initially 
increases and reaches a plateau for large systems. Both plots 
are in logarithmic scale. 



We have studied random (RP) and Achlioptas pro- 
cesses (APSR and APPR) for site percolation on 2D 
square lattices. We report the critical point for the APSR 
as pc = 0.695 and pc = 0.756 for APPR. We have shown 
that the order parameter Pmax changes more abruptly for 
the case of AP processes compared to the RP case. Using 
standard finite size scaling techniques, we find indications 
that "explosive" site and bond percolation may belong to 
different universality classes. A qualitatively similar be- 
havior was also observed for A /TV, where the exponents 
for "explosive" site and bond percolation are different. 
Finally, the quantity x/^ shows at Pc a power-law decay 
for the case of RP, while it appears to be nearly constant 



for large systems in the case of AP models (APSR and 
APPR). 

We also have studied the direct and two variants of 
a reverse AP process. We have shown that, although 
in finite size systems, the reverse processes are different 
from the direct process, they become both equivalent in 
the thermodynamic limit. 
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FIG. 7: (a) Plot of the direct (black squares) and reverse 
(red dots)APl for a 700 x 700 lattice. The hysteresis loop is 
evident, (b) The area of the hysteresis loop AH as a function 
of 1/L (logarithmic scale). Black squares are for reverse API 
and red dots are for reverse AP2, as described in Sec. II. For 
reverse API, AH ~ (l/L)^-^^^ and for reverse AP2 AH ~ 
(l/L)°-°^^ It is clear that the reverse and forward processes 
become equivalent in the thermodynamic limit {L oo). 
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